Ground extraction method for 3d point clouds of outdoor scenes based on gaussian process regression

ABSTRACT

A ground extraction method for 3D point clouds of outdoor scenes based on Gaussian process regression, including: (1) obtaining the 3D point cloud of an outdoor scene, (2) building the neighborhood of the 3D point cloud, (3) calculating the covariance matrices and normal vectors of the 3D point cloud, (4) classifying the 3D point cloud according to its neighborhood shape, (5) extracting the initial ground G s , (6) segmenting the initial ground, (7) 2D Gaussian process regression, (8) finding the neighborhood N LG     k    of each ground fragment LG k , and (9) extracting the final ground G e .

CROSS-REFERENCE TO RELATED APPLICATIONS

This application is a continuation-in-part of International Patent Application No PCT/CN2021/071055 with an international filing date of Jan. 11, 2021, designating the United States, now pending, and further claims priority benefits to Chinese Patent Application No. 202010510160.6, filed Jun. 8, 2020. The contents of all of these specifications are incorporated herein by reference.

BACKGROUND OF THE INVENTION Field of the Invention

The invention relates to a ground extraction method for 3D point clouds, more specifically, to a ground extraction method for 3D point clouds of outdoor scenes based on Gaussian process regression.

Description of the Related Art

With the development of 3D laser ranging technology, 3D point clouds are widely used in reverse engineering, industrial inspection, autonomous navigation, and so on. As the basis of the above applications, the 3D point cloud processing technology plays a vital role. In the 3D point cloud processing technology, feature extraction of 3D point clouds is an important technique, especially ground feature extraction of 3D point clouds of outdoor scenes, which plays an important role in segmentation and recognition of outdoor scenes, path planning of mobile robots, and so on.

In the field of autonomous navigation of mobile robots, the ground extraction of 3D point clouds of outdoor scenes is the premise of path planning of mobile robots. A complete ground provides passable areas for mobile robots, improves the passing capacity of mobile robots, and ensures the safety of mobile robots during their moving. In the field of outdoor scene analysis, since outdoor scenes are extremely complex, various objects are involved, such as buildings, trees, vehicles, people, and so on. To facilitate outdoor scene analysis, it is necessary to effectively segment 3D point clouds of outdoor scenes. Since the ground is the background of outdoor scenes, ground extraction contributes to separating the objects on the ground, which facilitates the object segmentation and scene analysis.

Now, the commonly used ground extraction method of 3D point clouds is the random sampling consistency (RANSAC) algorithm, which regards the ground as the largest surface in an outdoor scene. The method has a good performance for the flat and large ground. However, for some complex outdoor scenes, the ground is fragmentary and fluctuant. The method cannot guarantee the integrity and accuracy of ground extraction.

SUMMARY OF THE INVENTION

In view of the above-described problems, it is one objective of the invention to provide a ground extraction method for 3D point clouds of outdoor scenes based on Gaussian process regression. For an outdoor scene, the method first uses a laser rangefinder to obtain the 3D point cloud of an outdoor scene, which is essentially a set of discrete points in 3D space. Then, through the ground extraction method, the ground point cloud data are accurately and completely extracted from the 3D point cloud of the outdoor scene. The method solves the problem of incomplete and inaccurate ground extraction caused by the factors, such as complex outdoor scenes, fragmentary ground, and fluctuant ground, so as to improve the accuracy and integrity of ground extraction in outdoor scenes. The method has good performance on ground extraction.

To achieve the above objectives, in accordance with one embodiment of the invention, there is provided a ground extraction method of 3D point clouds of outdoor scenes based on Gaussian process regression, which includes the following steps.

Step 1: Obtaining the 3D point cloud of an outdoor scene. The 3D point cloud of an outdoor scene (a set of discrete points) is collected by using a laser rangefinder.

Step 2: Building the neighborhood of the 3D point cloud. The structure tree of the 3D point cloud is constructed by using the KD-tree algorithm. According to the coordinates of the discrete points in the 3D point cloud, the 3D point cloud is divided into different spatial regions. In the process of the neighborhood construction, the spatial address information is used to search the neighboring points. Then, the neighborhood N={p_(i)=(x_(i), y_(i), z_(i))|1≤i≤n_(n)} of a given point p=(x,y,z) in the 3D point cloud is built, where p_(i) is a neighboring point of the given point p=(x,y,z), i is the serial number of the neighboring points p_(i) and n_(n) is the number of the neighboring points p_(i).

Step 3: Calculating the covariance matrices and normal vectors of the 3D point cloud. For the given point p=(x,y,z) in the 3D point cloud, its covariance matrix M is constructed by using its neighborhood N={p_(i)=(x_(i),y_(i),z_(i))|1≤i≤n_(n)}. Then, the eigenvalues λ₁, λ₂, λ₃ and eigenvectors v₁,v₂,v₃ of the covariance matrix M and the normal vector n of the given point p are computed. This step includes the following substeps:

(a) By using the neighborhood relationship built in step 2, the neighborhood N={p_(i)=(x_(i),y_(i),z_(i))|1≤i≤n_(n)} of the given point p=(x,y,z) is built.

(b) The covariance matrix M of the neighborhood N of the given point p is computed as

M=Σ _(i=1) ^(n) ^(n) (p _(i) −p)(p _(i) −p)^(T)  (1),

where T is a vector transpose symbol that transposes a column vector to a row vector.

(c) The eigenvalues λ₁, λ₂, λ₃ (λ₁<λ₂<λ₃) and the corresponding eigenvectors v₁, v₂, v₃ of the covariance matrix M is computed.

(d) The normal vector n of the given point p is computed as the eigenvector v₁ corresponding to the minimum eigenvalue λ₁.

(e) For each point in the 3D point cloud, its eigenvalues, eigenvectors, and normal vector are computed by using the substeps (a)-(d).

Step 4: Classifying the 3D point cloud according to its neighborhood shape. The neighborhood shape of the given point p is determined by using the relationship among the eigenvalues λ₁, λ₂, λ₃ of the covariance matrix M of the given point p. The 3D point cloud is divided into three classes, the scattered points C_(p), the linear points C_(l), and the surface points C_(s). This step includes the following substeps:

(a) If λ₁≈λ₂≈λ₃, i.e. λ₃/λ₂≤8 and λ₂/λ₁≤8, the given point p and its neighboring points p_(i) present a scattered shape, and the given point p is categorized as a scattered point.

(b) If λ₁≈λ₂<<λ₃, i.e. λ₃/λ₂>8 and λ₂/λ₁≤8, the given point p and its neighboring points p_(i) present a linear shape, and the given point p is categorized as a linear point.

(c) If λ₁<<λ₂≈λ₃, i.e. λ₃/λ₂≤8 and λ₂/λ₁>8, the given point p and its neighboring points p_(i) present a surface shape, and the given point p is categorized as a surface point.

(d) For each point in the 3D point cloud, the substeps (a)-(c) are performed. Then, the 3D point cloud is divided into three classes, the scattered points C_(p), the linear points C_(l), and the surface points C_(s).

Step 5: Extracting the initial ground G_(s). For each surface point in the set C_(s), its normal vector is put on a unit ball S, which is called the normal ball. By using the mean shift algorithm, the normal vectors of the surface points C_(s) are clustered on the normal ball S Then, the surface points are divided into several surface regions F_(j). Finally, the initial ground G_(s) is extracted from the surface regions. This step includes the following substeps.

(a) For each surface point in C_(s), its normal vector is put on a unit ball S, which is called the normal ball. That is, the terminal point of the normal vector is located on the surface of the normal ball S, and the initial point of the normal vector is located at the center of the normal ball S.

(b) The mean shift algorithm is chosen to cluster the normal vectors of the surface points C_(s) on the normal ball S. Then, the normal vectors of the surface points C_(s) are divided into several groups. And, the surface points C_(s) are also divided into several surface regions F_(j) correspondingly, 1≤j≤m, where j is the serial number of the surface regions and m is the number of the surface regions.

(c) For each surface region, its average elevation h _(j) and average normal vector n _(j) are calculated. If the average elevation h _(j) and average normal vector n _(j) meet the conditions: h _(j)≤0.5 m and −5°≤α _(j)≤+5°, where α _(j) is the included angle between the average normal vector n _(j) and the vertical direction, the surface region F_(j) is considered as a part of the initial ground G_(s). This method is applied to each region F_(j), and then the initial ground G_(s) is obtained.

Step 6: Segmenting the initial ground. The initial ground point p_(t) is the point in the initial ground G_(s) By using the K-means algorithm, the initial ground G_(s) is divided into K ground fragments according to the distance between the initial ground point p_(t) and the center point p_(k) of each ground fragment. Each ground fragment is denoted by LG_(k)={p_(ks)=(x_(ks),y_(ks),z_(ks))|1≤s≤n_(k)}, 1≤k≤K, where p_(ks) is the point in the ground fragment LG_(k), k is the serial number of the ground fragments, K is the number of the ground fragments, s is the serial number of the points in the ground fragment LG_(k), and n_(k) is the number of the points in the ground fragment LG_(k). This step includes the following substeps:

(a) By using the farthest point sampling method, K points are determinded as the initial center points p_(k) of the ground fragments, 1≤k≤K. K is set at the beginning.

(b) The distance between each initial ground point p_(t) (1≤t≤n_(s)) and the center point p_(k) of each ground fragment is calculated, where t is the serial number of the points in the initial ground G_(s) and n_(s) is the number of the points in the initial ground G_(s). Each initial ground point p_(t) is assigned to its closest center point of the ground fragment.

(c) When all the initial ground points p_(t) are assigned, the center point p_(k) of each ground fragment is recalculated according to the points p_(ks) in the ground fragment LG_(k) by

$\begin{matrix} {{p_{k} = {\frac{1}{n_{k}}{\sum\limits_{s = 1}^{n_{k}}p_{ks}}}},} & (2) \end{matrix}$

where p_(k) is the new center point of the ground fragment.

(d) The substeps (b)-(c) are repeated until p_(k) does not change again. Finally, the points p_(ks) assigned to p_(k) are grouped as one ground fragment. The initial ground G_(s) is divided into K ground fragments LG_(k) (1≤k≤K).

Step 7: 2D Gaussian process regression. 2D Gaussian process regression is used to build the probabilistic model for each ground fragment LG_(k). The ground fragment LG_(k) is used as the training samples to train the 2D Gaussian process regression model. The hyperparameters l_(k), σ_(k) ² and σ_(n) ² of the model are solved by using the conjugate gradient method. This step includes the following substeps:

(a) A discrete function z_(ks)=ƒ_(d)(t_(ks)), t_(ks)=[x_(ks), y_(ks)]^(T) is defined on the ground fragment LG_(k). The value of the discrete function ƒ_(d)(t_(ks)) at a location t_(ks) represents the random variable in the Gaussian process.

(b) The Gaussian process is specified by its mean function m(t_(ks)) and covariance function k(t_(ks), t_(kt)), t_(ks), t_(kt)∈{t_(k1), t_(k2), . . . , t_(kn) _(k) } as

$\begin{matrix} \left\{ {\begin{matrix} {{m\left( t_{ks} \right)} = {\left( {f_{d}\left( t_{ks} \right)} \right)}} \\ {{k\left( {t_{ks},t_{kt}} \right)} = {\left\lbrack {\left( {{f_{d}\left( t_{ks} \right)} - {m\left( t_{ks} \right)}} \right)\left( {{f_{d}\left( t_{kt} \right)} - {m\left( t_{kt} \right)}} \right)} \right\rbrack}} \end{matrix},} \right. & (3) \end{matrix}$

Thus, the Gaussian process can be written as ƒ_(d)/(t_(ks))˜

(m(t_(ks)), k(t_(ks), t_(kt))).

(c) In consideration of the noise ε, the discrete function is changed to z_(ks)=ƒ_(d)(t_(ks))+ε. Assume additive identically distributed independent Gaussian noise ε with the variance σ_(n) ². The joint prior distribution of the observed values z_(k)=[z_(k1),z_(k2), . . . , z_(kn) _(k) ]^(T) is

z _(k) ˜N(0,K(T _(k) ,T _(k))+σ_(n) ² I)  (4),

where T_(k)=[t_(k1),t_(k2), . . . ,t_(kn) _(k) ] and K(T_(k),T_(k))=(k_(st)) is a n_(k)×n_(k) covariance matrix of the training inputs T_(k). The element k_(st) is used to measure the correlation between t_(ks) and t_(kt), and k_(st) is the squared exponential covariance function, which has the following form

$\begin{matrix} {{k_{st} = {\sigma_{k}^{2}\exp\left\{ {{- \frac{1}{2l_{k}^{2}}}{{t_{ks} - t_{kt}}}^{2}} \right\}}},} & (5) \end{matrix}$

where l_(k) is the length-scale, σ_(k) ² is the signal variance, and 1≤s, t≤n_(k)(s≠t).

(d) The joint prior distribution of the observed values z_(k) and the predicted value ƒ_(*) at the test location t_(*)=[x_(*),y_(*)]^(T) is given by

$\begin{matrix} {\begin{bmatrix} z_{k} \\ f_{*} \end{bmatrix} \sim {N\left( {0,\begin{bmatrix} {{K\left( {T_{k},T_{k}} \right)} + {\sigma_{n}^{2}I}} & {K\left( {T_{k},t_{*}} \right)} \\ {K\left( {t_{*},T_{k}} \right)} & {K\left( {t_{*},t_{*}} \right)} \end{bmatrix}} \right)}} & (6) \end{matrix}$

where K(T_(k),t_(*))=K(t_(*),T_(k))^(T) is the n_(k)×1 covariance matrix of the training inputs T_(k) and the test input t_(*), K(t_(*),t_(*)) is the covariance of the test input t_(*), and I is the n_(k)×n_(k) unit matrix. Then, the posterior distribution of the predicted value ƒ_(*), namely the 2D Gaussian process regression model of the ground fragment LG_(k), is given by

ƒ_(*) |T _(k) ,z _(k) ,t _(*) ˜N(m(ƒ_(*)),cov(ƒ_(*)))  (7),

where

m(ƒ_(*))=K(t _(*) ,T _(k))[K(T _(k) ,T _(k))+σ_(n) ² I]⁻¹ z _(k)  (8),

cov(ƒ_(*))=K(t _(*) ,t _(*))−K(t _(*) ,T _(k))[K(T _(k) ,T _(k))+σ_(n) ² I]⁻¹ K(T _(k) ,t _(*))  (9),

m(ƒ_(*)) and cov(ƒ_(*)) are the mean and variance of ƒ_(*) at the test location t_(*).

(e) The points in the ground fragment LG_(k) are used as the training samples. The negative logarithm likelihood function of the conditional probability of the training samples is established. And, the partial derivatives of the hyperparameters l_(k), σ_(k) ² and σ_(n) ² are calculated. Then, the conjugate gradient method is used to minimize the partial derivatives to obtain the optimal solution of the hyperparameters.

Step 8: Finding the neighborhood N_(LG) _(k) of each ground fragment LG_(k). For each ground fragment LG_(k), the distance between each point p_(ks) (p_(ks)∈LG_(k)) in the ground fragment LG_(k) and each point p_(a) in the 3D point cloud P={p_(a)=(x_(a), y_(a), z_(a))|1≤a≤n_(a)} of the outdoor scene is computed. If the distance is less than the threshold r_(k), p_(a) is regarded as a neighboring point of the ground fragment LG_(k), where a is the serial number of the points in the 3D point cloud P of the outdoor scene and n_(a) is the number of the points in the 3D point cloud P of the outdoor scene. This step includes the following substeps.

(a) The neighborhood N_(p) _(ks) ={p_(ks) ^(b)|1≤b≤n_(ks),p_(ks) ^(b)∈P,|p_(ks) ^(b)−p_(ks)|≤r_(k)} of each points p_(ks) in the ground fragment LG_(k) is computed, where n_(ks) is the number of the neighboring points of p_(ks), r_(k)=0.25√{square root over (n_(k))}d_(k), d_(k)=Σ_(s=1) ^(n) ^(k) |p_(ks)−p _(ks)|/n_(k), and p _(ks) is the closest point to p_(ks).

(b) All the neighborhoods N_(p) _(ks) compose the neighborhood of the ground fragment LG_(k). Let N_(LG) _(k) ={{circumflex over (p)}_(ks′)=({circumflex over (x)}_(ks′),ŷ_(ks′),{circumflex over (z)}_(ks′))|1≤s′≤{circumflex over (n)}_(k)}, where {circumflex over (p)}_(ks′) is the neighboring point of the ground fragment LG_(k), s′ is the serial number of the neighboring points, and {circumflex over (n)}_(k) is the number of the neighboring points.

Step 9: Extracting the final ground G_(e). The neighborhood of each ground fragment LG_(k) is substituted into the 2D Gaussian process regression model of the ground fragment LG_(k). For each point {circumflex over (p)}_(ks′)∈N_(LG) _(k) , the predicted mean m({circumflex over (ƒ)}_(ks′)) and variance cov({circumflex over (ƒ)}_(ks′)) at the test location {circumflex over (t)}_(ks′)=[{circumflex over (x)}_(ks′),ŷ_(ks′)]^(T) are calculated. If they meet the conditions: cov({circumflex over (ƒ)}_(ks′))≤σ ² and |m({circumflex over (ƒ)}_(ks′))−{circumflex over (z)}_(ks′)|/√{square root over (σ_(k) ²+σ_(n) ²)}≤d, where σ ²=0.05 is the threshold of the variance and d=0.2 is the threshold of the Mahalanobis distance, the point {circumflex over (p)}_(ks′) is regarded as a ground point. By using this method for each ground fragment LG_(k) to check the points in its neighborhood, the ground points in the neighborhood N_(LG) _(k) of each ground fragment LG_(k) is found. Then, the final ground G_(e) is obtained.

The beneficial effect of the invention is a ground extraction method for 3D point clouds of outdoor scenes based on Gaussian process regression, which includes the following steps. (1) obtaining the 3D point cloud of an outdoor scene, (2) building the neighborhood of the 3D point cloud, (3) calculating the covariance matrices and normal vectors of the 3D point cloud, (4) classifying the 3D point cloud according to its neighborhood shape, (5) extracting the initial ground G_(s), (6) segmenting the initial ground, (7) 2D Gaussian process regression, (8) finding the neighborhood N_(LG) _(k) of each ground fragment LG_(k), and (9) extracting the final ground G_(e). Compared with the existing methods, the invention adopts the idea of layer-by-layer extraction First, through the eigenvalue analysis of covariance matrices, the whole surface area (the set of surface points C_(s)) in the 3D point cloud of the outdoor scene is extracted. Then, by building the normal ball S of the set of surface points C_(s) and clustering the normal vectors on the normal ball S, the whole surface area (the set of surface points C_(s)) is divided into several surface regions F_(j). Through the combination of normal vectors and elevations, the initial ground G_(s) is obtained from surface regions F_(j). The K-means clustering algorithm is used to segment the initial ground G_(s) into several more compact ground fragments LG_(k). Finally, the final ground G_(e) is obtained by Gaussian process regression of the ground fragments LG_(k). The idea of layer-by-layer extraction can make the ground extraction of 3D point clouds more complete and accurate. Especially, when the outdoor scene very is complex, the initial ground is divided into several compact ground fragments, the Gaussian process regression is then carried out respectively, and the final ground is obtained by OR operation at last. This contributes to extracting the fragmentary and fluctuant ground. Therefore, the method of the invention effectively solves the problem of incomplete and inaccurate ground extraction caused by the factors, such as complex outdoor scenes, fragmentary ground, and fluctuant ground. The method has good performance on ground extraction.

BRIEF DESCRIPTION OF THE DRAWINGS

The invention is described hereinbelow with reference to the accompanying drawings, in which:

FIG. 1 is the flow chart of the method of the invention.

FIG. 2 is the 3D point cloud of an outdoor scene.

FIG. 3 is the result of the surface point extraction.

FIG. 4 is the result of the normal ball construction.

FIG. 5 is the result of the normal vector clustering.

FIG. 6 is the result of the initial ground extraction.

FIG. 7 is the result of the initial ground segmentation.

FIG. 8 is a ground fragment and its neighborhood.

FIG. 9 is the result of the final ground extraction.

DETAILED DESCRIPTION OF THE DRAWINGS

To further illustrate the invention, embodiments detailing a ground extraction method for 3D point clouds are described below. It should be noted that the following embodiments are intended to describe and not to limit the invention.

The invention will be further explained with the figures.

As shown in FIG. 1, a ground extraction method for 3D point clouds of outdoor scenes based on Gaussian process regression. Its characteristics include the following steps:

Step 1: Obtaining the 3D point cloud of an outdoor scene. The 3D point cloud of an outdoor scene (a set of discrete points) is collected by using a laser rangefinder. As shown in FIG. 2, the whole outdoor scene consists of about 100000 points, including the ground, trees, grass, buildings, vehicles, people, etc.

Step 2: Building the neighborhood of the 3D point cloud. The structure tree of the 3D point cloud is constructed by using the KD-tree algorithm. According to the coordinates of the discrete points in the 3D point cloud, the 3D point cloud is divided into different spatial regions. In the process of the neighborhood construction, the spatial address information is used to search the neighboring points. Then, the neighborhood N={p_(i)=(x_(i), y_(i), z_(i))|1≤i≤n_(n)} of a given point p=(x,y,z) in the 3D point cloud is built, where p_(i) is a neighboring point of the given point p=(x,y,z), i is the serial number of the neighboring points p_(i) and n_(n) is the number of the neighboring points p_(i).

Step 3: Calculating the covariance matrices and normal vectors of the 3D point cloud. For the given point p=(x,y,z) in the 3D point cloud, its covariance matrix M is constructed by using its neighborhood N={p_(i)=(x_(i), y_(i),z_(i))|1≤i≤n_(n)}. Then, the eigenvalues Δ₁,λ₂,λ₃ and eigenvectors v₁,v₂,v₃ of the covariance matrix M and the normal vector n of the given point p are computed. This step includes the following substeps:

(a) By using the neighborhood relationship built in step 2, the neighborhood N={p_(i)=(x_(i),y_(i),z_(i))|1≤i≤n_(n)} of the given point p=(x,y,z) is built.

(b) The covariance matrix M of the neighborhood N of the given point p is computed as

M=Σ _(i=1) ^(n) ^(n) (p _(i) −p)(p _(i) −p)^(T)  (1),

where T is a vector transpose symbol that transposes a column vector to a row vector.

(c) The eigenvalues λ₁, λ₂, λ₃ (λ₁<λ₂<λ₃) and the corresponding eigenvectors v₁, v₂, v₃ of the covariance matrix M is computed.

(d) The normal vector n of the given point p is computed as the eigenvector v₁ corresponding to the minimum eigenvalue λ₁.

(e) For each point in the 3D point cloud, its eigenvalues, eigenvectors, and normal vector are computed by using the substeps (a)-(d).

Step 4: Classifying the 3D point cloud according to its neighborhood shape. The neighborhood shape of the given point p is determined by using the relationship among the eigenvalues λ₁, λ₂, λ₃ of the covariance matrix M of the given point p. The 3D point cloud is divided into three classes: the scattered points C_(p), the linear points C_(l), and the surface points C_(s). This step includes the following substeps:

(a) If λ₁≈λ₂≈λ₃, i.e. λ₃/λ₂≤8 and λ₂/λ₁≤8, the given point p and its neighboring points p_(i) present a scattered shape, and the given point p is categorized as a scattered point.

(b) If λ₁≈λ₂<<λ₃, i.e. λ₃/λ₂>8 and λ₂/λ₁≤8, the given point p and its neighboring points p_(i) present a linear shape, and the given point p is categorized as a linear point.

(c) If λ₁<<λ₂≈λ₃, i.e. λ₃/λ₂≤8 and λ₂/λ₁>8, the given point p and its neighboring points p_(i) present a surface shape, and the given point p is categorized as a surface point.

(d) For each point in the 3D point cloud, the substeps (a)-(c) are performed. Then, the 3D point cloud is divided into three classes, the scattered points C_(p), the linear points C_(l), and the surface points C_(s). The extraction result of the surface points C_(s) of the outdoor scene is shown in FIG. 3.

Step 5: Extracting the initial ground G_(s) For each surface point in the set C_(s), its normal vector is put on a unit ball S, which is called the normal ball. By using the mean shift algorithm, the normal vectors of the surface points C_(s) are clustered on the normal ball S. Then, the surface points are divided into several surface regions F_(j). Finally, the initial ground G_(s) is extracted from the surface regions. This step includes the following substeps:

(a) For each surface point in C_(s), its normal vector is put on a unit ball S, which is called the normal ball. As shown in FIG. 4, the terminal point of the normal vector is located on the surface of the normal ball S, and the initial point of the normal vector is located at the center of the normal ball S.

(b) The mean shift algorithm is chosen to cluster the normal vectors of the surface points C_(s) on the normal ball S. As shown in FIG. 5, the normal vectors of the surface points C_(s) are divided into several groups. And, the surface points C_(s) are also divided into several surface regions F_(j) correspondingly, 1≤j≤m, where j is the serial number of the surface regions and m is the number of the surface regions.

(c) For each surface region, its average elevation h _(j) and average normal vector n _(j) are calculated. If the average elevation h _(j) and average normal vector n _(j) meet the conditions. h _(j)≤0.5m and −5°≤α _(j)≤+5°, where α _(j) is the included angle between the average normal vector n _(j) and the vertical direction, the surface region F_(j) is considered as a part of the initial ground G_(s). As shown in FIG. 6, this method is applied to each region F_(j), and then the initial ground G_(s) is obtained.

Step 6: Segmenting the initial ground. The initial ground point p_(t) is the point in the initial ground G_(s). By using the K-means algorithm, the initial ground G_(s) is divided into K ground fragments according to the distance between the initial ground point p_(t) and the center point p_(k) of each ground fragment. Each ground fragment is denoted by LG_(k)={p_(ks)=(x_(ks),y_(ks),z_(ks))|1≤s≤n_(k)}, 1≤k≤K, where p_(ks) is the point in the ground fragment LG_(k), k is the serial number of the ground fragments, K is the number of the ground fragments, s is the serial number of the points in the ground fragment LG_(k), and n_(k) is the number of the points in the ground fragment LG_(k). This step includes the following substeps:

(a) By using the farthest point sampling method, K points are determined as the initial center points p_(k) of the ground fragments, 1≤k≤K. K is set at the beginning.

(b) The distance between each initial ground point p_(t) (1≤t≤n_(s)) and the center point p_(k) of each ground fragment is calculated, where t is the serial number of the points in the initial ground G_(s) and n_(s) is the number of the points in the initial ground G_(s). Each initial ground point p_(t) is assigned to its closest center point of the ground fragment.

(c) When all the initial ground points p_(t) are assigned, the center point p_(k) of each ground fragment is recalculated according to the points p_(ks) in the ground fragment LG_(k) by

$\begin{matrix} {{p_{k} = {\frac{1}{n_{k}}{\sum\limits_{s = 1}^{n_{k}}p_{ks}}}},} & (2) \end{matrix}$

where p_(k) is the new center point of the ground fragment.

(d) The substeps (b)-(c) are repeated until p_(k) does not change again. Finally, the points p_(k) assigned to p_(k) are grouped as one ground fragment. The initial ground G_(s) is divided into K ground fragments LG_(k)(1≤k≤K), as shown in FIG. 7.

Step 7: 2D Gaussian process regression. 2D Gaussian process regression is used to build the probabilistic model for each ground fragment LG_(k). The ground fragment LG_(k) is used as the training samples to train the 2D Gaussian process regression model. The hyperparameters l_(k), σ_(k) ² and σ_(n) ² of the model are solved by using the conjugate gradient method. This step includes the following substeps:

(a) A discrete function z_(ks)=ƒ_(d)(t_(ks)), t_(ks)=[x_(ks), y_(ks)]^(T) is defined on the ground fragment LG_(k). The value of the discrete function ƒ_(d)(t_(ks)) at a location t_(ks) represents the random variable in the Gaussian process.

(b) The Gaussian process is specified by its mean function m(t_(ks)) and covariance function k(t_(ks), t_(kt)), t_(ks), t_(kt)∈{t_(k1), t_(k2), . . . , t_(kn) _(k) } as

$\begin{matrix} \left\{ {\begin{matrix} {{m\left( t_{ks} \right)} = {\left( {f_{d}\left( t_{ks} \right)} \right)}} \\ {{k\left( {t_{ks},t_{kt}} \right)} = {\left\lbrack {\left( {{f_{d}\left( t_{ks} \right)} - {m\left( t_{ks} \right)}} \right)\left( {{f_{d}\left( t_{kt} \right)} - {m\left( t_{kt} \right)}} \right)} \right\rbrack}} \end{matrix}.} \right. & (3) \end{matrix}$

Thus, the Gaussian process can be written as

ƒ_(d)(t _(ks))˜

(m(t _(ks)),k(t _(ks) ,t _(kt))).

(c) In consideration of the noise ε, the discrete function is changed to z_(ks)=ƒ_(d)(t_(ks))+ε. Assume additive identically distributed independent Gaussian noise ε with the variance σ_(n) ². The joint prior distribution of the observed values z_(k)=[z_(k1), z_(k2), . . . , z_(kn) _(k) ]^(T) is

z _(k) ˜N(0,K(T _(k) ,T _(k))+σ_(n) ² I)  (4),

where T_(k)=[t_(k1), t_(k2), . . . , t_(kn) _(k) ] and K(T_(k),T_(k))=(k_(st)) is a n_(k)×n_(k) covariance matrix of the training inputs T_(k). The element k_(st) is used to measure the correlation between t_(ks) and t_(kt), and k_(st) is the squared exponential covariance function, which has the following form

$\begin{matrix} {{k_{st} = {\sigma_{k}^{2}\exp\left\{ {{- \frac{1}{2l_{k}^{2}}}{{t_{ks} - t_{kt}}}^{2}} \right\}}},} & (5) \end{matrix}$

where l_(k) is the length-scale, σ_(k) ² of is the signal variance, and 1≤s, t≤n_(k)(s≠t).

(d) The joint prior distribution of the observed values z_(k) and the predicted value ƒ_(*) at the test location t_(*)=[x_(*),y_(*)]^(T) is given by

$\begin{matrix} {{\begin{bmatrix} z_{k} \\ f_{*} \end{bmatrix} \sim {N\left( {0,\begin{bmatrix} {{K\left( {T_{k},T_{k}} \right)} + {\sigma_{n}^{2}I}} & {K\left( {T_{k},t_{*}} \right)} \\ {K\left( {t_{*},T_{k}} \right)} & {K\left( {t_{*},t_{*}} \right)} \end{bmatrix}} \right)}},} & (6) \end{matrix}$

where K(T_(k),t_(*))=K(t_(*),T_(k))^(T) is the n_(k)×1 covariance matrix of the training inputs T_(k) and the test input t_(*), K(t_(*),t_(*)) is the covariance of the test input t_(*), and I is the n_(k)×n_(k) unit matrix. Then, the posterior distribution of the predicted value ƒ_(*), namely the 2D Gaussian process regression model of the ground fragment LG_(k), is given by

ƒ_(*) |T _(k) ,z _(k) ,t _(*) ˜N(m(ƒ_(*)),cov(ƒ_(*)))  (7),

where

m(ƒ_(*))=K(t _(*) ,T _(k))[K(T _(k) ,T _(k))+σ_(n) ² I]⁻¹ z _(k)  (8),

cov(ƒ_(*))=K(t _(*) ,t _(*))−K(t _(*) ,T _(k))[K(T _(k) ,T _(k))+σ_(n) ² I]⁻¹ K(T _(k) ,t _(*))  (9),

m(ƒ_(*)) and cov(ƒ_(*)) are the mean and variance of ƒ_(*) at the test location t_(*).

(e) The points in the ground fragment LG_(k) are used as the training samples. The negative logarithm likelihood function of the conditional probability of the training samples is established And, the partial derivatives of the hyperparameters l_(k), σ_(k) ² and σ_(n) ² are calculated. Then, the conjugate gradient method is used to minimize the partial derivatives to obtain the optimal solution of the hyperparameters.

Step 8: Finding the neighborhood N_(LG) _(k) of each ground fragment LG_(k). For each ground fragment LG_(k), the distance between each point p_(ks) (p_(ks) ∈LG_(k)) in the ground fragment LG_(k) and each point p_(a) in the 3D point cloud P={p_(a)=(x_(a), y_(a), z_(a))|1≤a≤n_(a)} of the outdoor scene is computed. If the distance is less than the threshold r_(k), p_(a) is regarded as a neighboring point of the ground fragment LG_(k), where a is the serial number of the points in the 3D point cloud P of the outdoor scene and n_(a) is the number of the points in the 3D point cloud P of the outdoor scene. This step includes the following substeps.

(a) The neighborhood N_(p) _(ks) ={p_(ks) ^(b)|1≤b≤n_(ks), p_(ks) ^(b)∈P,|p_(ks) ^(b)−p_(ks)|≤r_(k)} of each points p_(ks) in the ground fragment LG_(k) is computed, where n_(ks) is the number of the neighboring points of p_(ks), r_(k)=0.25√{square root over (n_(k))}d_(k), d_(k)=Σ_(s=1) ^(n) ^(k) |p_(ks)−p _(ks)|/n_(k), and p _(ks) is the closest point to p_(ks).

(b) All the neighborhoods N_(p) _(ks) compose the neighborhood N_(LG) _(k) of the ground fragment LG_(k). Let N_(LG) _(k) ={{circumflex over (p)}_(ks′)=({circumflex over (x)}_(ks′),ŷ_(ks′),{circumflex over (z)}_(ks′))|1≤s′≤{circumflex over (n)}_(k)}, where {circumflex over (p)}_(ks′) is the neighboring point of the ground fragment LG_(k), s′ is the serial number of the neighboring points, and {circumflex over (n)}_(k) is the number of the neighboring points. The neighborhood N_(LG) _(k) of the ground fragment LG_(k) is shown in FIG. 8.

Step 9: Extracting the final ground G_(e) The neighborhood N_(LG) _(k) of each ground fragment LG_(k) is substituted into the 2D Gaussian process regression model of the ground fragment LG_(k). For each point {circumflex over (p)}_(ks′)∈N_(LG) _(k) , the predicted mean m({circumflex over (ƒ)}_(ks′)) and variance cov({circumflex over (ƒ)}_(ks′)) at the test location {circumflex over (t)}_(ks′)=[{circumflex over (x)}_(ks′),ŷ_(ks′)]^(T) are calculated. If they meet the conditions: cov({circumflex over (ƒ)}_(ks′))≤σ ² and |m({circumflex over (ƒ)}_(ks′))−{circumflex over (z)}_(ks′)|/√{square root over (σ_(k) ²+σ_(n) ²)}≤d, where σ ²=0.05 is the threshold of the variance and 3=0.2 is the threshold of the Mahalanobis distance, the point {circumflex over (p)}_(ks′) is regarded as a ground point. By using this method for each ground fragment LG_(k) to check the points in its neighborhood, the ground points in the neighborhood N_(LG) _(k) of each ground fragment LG_(k) is found. Then, the final ground G_(e) is obtained as shown in FIG. 9.

Advantages of the ground extraction method for 3D point clouds of outdoor scenes based on Gaussian process regression according to the invention are described as follows. Without any function assumption, such as linear function and quadratic function, Gaussian process regression can well approximate the unknown function by considering the correlation between the function value and the given observation value, even in the case of the abrupt change. This feature facilitates the flexible handling of the irregular ground. The invention uses the idea of layer-by-layer extraction and Gaussian process regression to accurately and completely extract the ground from 3D point clouds of outdoor scenes, which effectively solves the problems of incomplete and inaccurate ground extraction caused by the factors, such as complex outdoor scene, fragmentary ground, and fluctuant ground. The method has good performance on ground extraction.

While particular embodiments of the invention have been shown and described, it will be obvious to those skilled in the art that changes and modifications may be made without departing from the invention in its broader aspects, and therefore, the aim in the appended claims is to cover all such changes and modifications as fall within the true spirit and scope of the invention 

What is claimed is:
 1. A ground extraction method for 3D point clouds of outdoor scenes based on Gaussian process regression, comprising the following steps: step 1: obtaining a 3D point cloud of an outdoor scene by using a laser rangefinder, wherein the 3D point cloud of the outdoor scene is a set of discrete points; step 2: building a neighborhood of the 3D point cloud, comprising: constructing a structure tree of the 3D point cloud by using a KD-tree algorithm, dividing the 3D point cloud into different spatial regions according to coordinates of the discrete points in the 3D point cloud, using spatial address information to search neighboring points in the process of the neighborhood construction, and building the neighborhood N={p_(i)=(x_(i),y_(i),z_(i))|1≤i≤n_(n)} of a given point p=(x,y,z) in the 3D point cloud, wherein p_(i) is a neighboring point of the given point p=(x,y,z), i is the serial number of the neighboring points p_(i), and n_(n) is the number of the neighboring points p_(i); step 3: calculating covariance matrices and normal vectors of the 3D point cloud, wherein for the given point p=(x,y,z) in the 3D point cloud, a covariance matrix M is constructed by using its neighborhood N={p_(i)=(x_(i),y_(i),z_(i))|1≤i≤n_(n)}, and eigenvalues λ₁,λ₂,λ₃ and eigenvectors v₁,v₂,v₃ of the covariance matrix M and a normal vector n of the given point p are computed; and step 3 comprises the following substeps: (a) building the neighborhood N={p_(i)=(x_(i),y_(i),z_(i))|1≤i≤n_(n)} of the given point p=(x,y,z) by using the neighborhood relationship built in step 2; (b) computing the covariance matrix M of the neighborhood N of the given point p as M=Σ _(i=1) ^(n) ^(n) (p _(i) −p)(p _(i) −p)^(T)  (1), wherein T is a vector transpose symbol that transposes a column vector to a row vector; (c) computing the eigenvalues λ₁, λ₂, λ₃ (λ₁<λ₂<λ₃) and the corresponding eigenvectors v₁, v₂, v₃ of the covariance matrix M, (d) computing the normal vector n of the given point p as the eigenvector v₁ corresponding to the minimum eigenvalue λ₁; and (e) for each point in the 3D point cloud, repeating the substeps (a)-(d), whereby computing eigenvalues, eigenvectors, and normal vector for each point in the 3D point cloud; step 4: classifying the 3D point cloud according to its neighborhood shape; wherein, a neighborhood shape of the given point p is determined by using the relationship among the eigenvalues λ₁, λ₂, λ₃ of the covariance matrix M of the given point p; the 3D point cloud is divided into three classes: the scattered points C_(p), the linear points C_(l), and the surface points C_(s); and step 4) comprises the following substeps: (a) categorizing the given point p as a scattered point if λ₁≈λ₂≈λ₃, i.e. λ₃/λ₂≤8 and λ₂/λ₁≤8, and the given point p and its neighboring points p_(i) present a scattered shape; (b) categorizing the given point p as a linear point if λ₁≈λ₂<<λ₃, i.e. λ₃/λ₂>8 and λ₂/λ₁≤8, and the given point p and its neighboring points p_(i) present a linear shape, (c) categorizing the given point pas a surface point if λ₁<<λ₂≈λ₃, i.e. λ₃/λ₂≤8 and λ₂/λ₁>8, and the given point p and its neighboring points p_(i) present a surface shape; and (d) repeating the substeps (a)-(c) for each point in the 3D point cloud, whereby dividing the 3D point cloud into three classes: the scattered points C_(p), the linear points C_(l), and the surface points C_(s); step 5: extracting an initial ground G_(s); wherein, for each surface point in the set C_(s), a normal vector thereof is put on a unit ball S, which is called a normal ball; by using a mean shift algorithm, the normal vectors of the surface points C_(s) are clustered on the normal ball S; then, the surface points are divided into several surface regions F_(j); finally, the initial ground G_(s) is extracted from the surface regions; and step 5 comprises the following substeps; (a) for each surface point in C_(s), putting its normal vector on a unit ball S, which is called the normal ball, wherein a terminal point of the normal vector is located on the surface of the normal ball S, and an initial point of the normal vector is located at the center of the normal ball S; (b) choosing the mean shift algorithm to cluster the normal vectors of the surface points C_(s) on the normal ball S; dividing the normal vectors of the surface points C_(s) into several groups; and dividing the surface points C_(s) into several surface regions F_(j) correspondingly, 1≤j≤m, where j is the serial number of the surface regions and m is the number of the surface regions; and (c) for each surface region, calculating its average elevation h _(j) and average normal vector n _(j); considering the surface region F_(j) as a part of the initial ground G_(s), if the average elevation h _(j) and average normal vector n _(j) meet the conditions: h _(j)≤0.5m and −5°≤α _(j)≤+5°, where α _(j) is the included angle between the average normal vector n _(j) and the vertical direction; and applying this method to each region F_(j) whereby obtaining the initial ground G_(s); step 6; segmenting the initial ground; wherein, the initial ground point p_(t) is the point in the initial ground G_(s); by using the K-means algorithm, the initial ground G_(s) is divided into K ground fragments according to the distance between the initial ground point p_(t) and the center point p_(k) of each ground fragment; each ground fragment is denoted by LG_(k)={p_(ks)=(x_(ks),y_(ks),z_(ks))|1≤s≤n_(k)}·1≤k≤K, where p_(ks) is the point in the ground fragment LG_(k), k is the serial number of the ground fragments, K is the number of the ground fragments, s is the serial number of the points in the ground fragment LG_(k), and n_(k) is the number of the points in the ground fragment LG_(k); and step 6 comprises the following substeps; (a) determining K points asinitial center points p_(k) of the ground fragments by using the farthest point sampling method, wherein 1≤k≤K, and K is set at the beginning; (b) calculating a distance between each initial ground point p_(t) (1≤t≤n_(s)) and the center point p_(k) of each ground fragment, where t is the serial number of the points in the initial ground G_(s) and n_(s) is the number of the points in the initial ground G_(s), and assigning each initial ground point p_(t) to its closest center point of the ground fragment; (c) when all the initial ground points p_(t) are assigned, recalculating the center point p_(k) of each ground fragment according to the points p_(ks) in the ground fragment LG_(k) by $\begin{matrix} {{p_{k} = {\frac{1}{n_{k}}{\sum\limits_{s = 1}^{n_{k}}p_{ks}}}},} & (2) \end{matrix}$ where p_(k) is the new center point of the ground fragment, and (d) repeating the substeps (b)-(c) until p_(k) does not change again; finally, grouping the points p_(ks) assigned to p_(k) as one ground fragment, and dividing the initial ground G_(s) into K ground fragments LG_(k) (1≤k≤K); step 7; performing 2D Gaussian process regression, wherein the 2D Gaussian process regression is used to build the probabilistic model for each ground fragment LG_(k); the ground fragment LG_(k) is used as the training samples to train the 2D Gaussian process regression model; the hyperparameters l_(k), σ_(k) ² and σ_(n) ² of the model are solved by using the conjugate gradient method; and step 7 comprises the following substeps: (a) defining a discrete function z_(ks)=ƒ_(d)(t_(ks)), t_(ks)=[x_(ks),y_(ks)]^(T) on the ground fragment LG_(k), wherein the value of the discrete function ƒ_(d)(t_(ks)) at a location t_(ks) represents the random variable in the Gaussian process; (b) specifying the Gaussian process by its mean function m(t_(ks)) and covariance function k(t_(ks), t_(kt)), t_(ks), t_(kt)∈{t_(k1), t_(k2), . . . , t_(kn) _(k) } as $\begin{matrix} \left\{ {\begin{matrix} {{m\left( t_{ks} \right)} = {\left( {f_{d}\left( t_{ks} \right)} \right)}} \\ {{k\left( {t_{ks},t_{kt}} \right)} = {\left\lbrack {\left( {{f_{d}\left( t_{ks} \right)} - {m\left( t_{ks} \right)}} \right)\left( {{f_{d}\left( t_{kt} \right)} - {m\left( t_{kt} \right)}} \right)} \right\rbrack}} \end{matrix}.} \right. & (3) \end{matrix}$  and writing Gaussian process as ƒ_(d)(t_(ks))˜

(m(t_(ks)), k(t_(ks), t_(kt))); (c) in consideration of the noise ε, changing the discrete function is to z_(ks)=ƒ_(d)(t_(ks))+ε; assuming additive identically distributed independent Gaussian noise ε with the variance σ_(n) ²; and obtaining the joint prior distribution of the observed values z_(k)=[z_(k1), z_(k2), . . . , z_(kn) _(k) ]^(T) as z _(k) ˜N(0,K(T _(k) ,T _(k))+σ_(n) ² I)  (4), where T_(k)=[t_(k1), t_(k2), . . . , t_(kn) _(k) ] and K(T_(k),T_(k))=(k_(st)) is a n_(k)×n_(k) covariance matrix of the training inputs T_(k); the element k_(st) is used to measure the correlation between t_(ks) and t_(kl), and k_(st) is the squared exponential covariance function, which has the following form $\begin{matrix} {{k_{st} = {\sigma_{k}^{2}\exp\left\{ {{- \frac{1}{2l_{k}^{2}}}{{t_{ks} - t_{kt}}}^{2}} \right\}}},} & (5) \end{matrix}$ where l_(k) is the length-scale, σ_(k) ² is the signal variance, and 1≤s, t≤n_(k)(s≠t); (d) giving the joint prior distribution of the observed values z_(k) and the predicted value ƒ_(*) at the test location t_(*)=[x_(*),y_(*)]^(T) by $\begin{matrix} {{\begin{bmatrix} z_{k} \\ f_{*} \end{bmatrix} \sim {N\left( {0,\begin{bmatrix} {{K\left( {T_{k},T_{k}} \right)} + {\sigma_{n}^{2}I}} & {K\left( {T_{k},t_{*}} \right)} \\ {K\left( {t_{*},T_{k}} \right)} & {K\left( {t_{*},t_{*}} \right)} \end{bmatrix}} \right)}},} & (6) \end{matrix}$ where K(T_(k),t_(*))=K(t_(*),T_(k))^(T) is the n_(k)×1 covariance matrix of the training inputs T_(k) and the test input t_(*), K(t_(*),t_(*)) is the covariance of the test input t_(*), and I is the n_(k)×n_(k) unit matrix; and giving the posterior distribution of the predicted value ƒ_(*), namely the 2D Gaussian process regression model of the ground fragment LG_(k), by ƒ_(*) |T _(k) ,z _(k) ,t _(*) ˜N(m(ƒ_(*)),cov(ƒ_(*)))  (7), where m(ƒ_(*))=K(t _(*) ,T _(k))[K(T _(k) ,T _(k))+σ_(n) ² I]⁻¹ z _(k)  (8), cov(ƒ_(*))=K(t _(*) ,t _(*))−K(t _(*) ,T _(k))[K(T _(k) ,T _(k))+σ_(n) ² I]⁻¹ K(T _(k) ,t _(*))  (9), m(ƒ_(*)) and cov(ƒ_(*)) are the mean and variance of ƒ_(*) at the test location t_(*); and (e) using the points in the ground fragment LG_(k) as the training samples, establishing the negative logarithm likelihood function of the conditional probability of the training samples; calculating the partial derivatives of the hyperparameters l_(k), σ_(k) ² and σ_(n) ²; and using the conjugate gradient method to minimize the partial derivatives to obtain the optimal solution of the hyperparameters; step 8: finding the neighborhood of each ground fragment LG_(k); wherein for each ground fragment LG_(k), a distance between each point p_(ks) (p_(ks)∈LG_(k)) in the ground fragment LG_(k) and each point p_(a) in the 3D point cloud P={p_(a)=(x_(a), y_(a), z_(a))|1≤a≤n_(a)} of the outdoor scene is computed; if the distance is less than the threshold r_(k), p_(a) is regarded as a neighboring point of the ground fragment LG_(k), where a is the serial number of the points in the 3D point cloud P of the outdoor scene and n_(a) is the number of the points in the 3D point cloud P of the outdoor scene, and step 8 comprises the following substeps: (a) computing the neighborhood N_(p) _(ks) ={p_(ks) ^(b)|1≤b≤n_(ks), p_(ks) ^(b)∈P, |p_(ks) ^(b)−p_(ks)|≤r_(k)} of each points p_(ks) in the ground fragment LG_(k), where n_(ks) is the number of the neighboring points of p_(ks), r_(k)=0.25√{square root over (n_(k))}d_(k), d_(k)=Σ_(s=1) ^(n) ^(k) |p_(ks)−p _(ks)|/n_(k), and p _(ks) is the closest point to p_(ks); and (b) composing, by all the neighborhoods N_(p) _(ks) , the neighborhood N_(LG) _(k) of the ground fragment LG_(k); and letting N_(LG) _(k) ={{circumflex over (p)}_(ks′)=({circumflex over (x)}_(ks′),ŷ_(ks′),{circumflex over (z)}_(ks′))|1≤s′≤{circumflex over (n)}_(k)}, where {circumflex over (p)}_(ks′) is the neighboring point of the ground fragment LG_(k), s′ is the serial number of the neighboring points, and {circumflex over (n)}_(k) is the number of the neighboring points; and step 9: extracting the final ground G_(e); wherein the neighborhood N_(LG) _(k) of each ground fragment LG_(k) is substituted into the 2D Gaussian process regression model of the ground fragment LG_(k); for each point {circumflex over (p)}_(ks′)∈N_(LG) _(k) , the predicted mean m({circumflex over (ƒ)}_(ks′)) and variance cov({circumflex over (ƒ)}_(ks′)) at the test location {circumflex over (t)}_(ks′)=[{circumflex over (x)}_(ks′),ŷ_(ks′)]^(T) are calculated, if they meet the conditions: cov({circumflex over (ƒ)}_(ks′))≤σ ² and |m({circumflex over (ƒ)}_(ks′))−{circumflex over (z)}_(ks′)|/√{square root over (σ_(k) ²+σ_(n) ²)}≤d, where σ ²=0.05 is the threshold of the variance and d=0.2 is the threshold of the Mahalanobis distance, the point {circumflex over (p)}_(ks′) is regarded as a ground point; by using this method for each ground fragment LG_(k) to check the points in its neighborhood, the ground points in the neighborhood N_(LG) _(k) of each ground fragment LG_(k) is found; then, the final ground G_(e) is obtained. 